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ABSTRACT 

This article introduces a family of analytical functions of the form x"K v (x), 
where K v is the incomplete Bessel function of the third kind. This family of 
functions can describe the density profile, projected and integrated light profiles 
and the gravitational potentials of galaxies. For the proper choice of parameters, 
these functions accurately approximate Sersic functions over a range of indices 
and are good fits to galaxy light profiles. With an additional parameter corre- 
sponding to a galaxy core radius, these functions can fit galaxy like M87 over 
a factor of 10 5 in radius. Unlike Sersic profiles, these functions have simple an- 
alytical 2-dimensional and 3-dimensional Fourier transforms, so they are easily 
convolved with spatially varying point spread function and are well suited for 
photometric and lensing analysis. We use these functions to estimate the effects 
of seeing on lensing measurements and show that high S/N measurements, even 
when the PSF is larger than the galaxy effective radius, should be able to recover 
accurate estimates of lensing distortions by weighting light in the outer isophotes 
that are less effected by seeing. 

1. Introduction 

Weak lensing observations have the potential to provide powerful new insights into the 



nature of dark energy and dark matter (see e.g., (Hoekstra & Jain 2008)) as well as directly 



the relationship between luminous and dark matter. 

Over the next few years, astronomers can anticipate very large, high quality photometric 
data. It is essential to develop image analysis techniques that can exploit this high quality 
data. The analysis techniques must be rapid and unbiased. Ideally, they should be nearly 



optimal and use most of the information in an astronomical image. While there has 
been significant progress in the past few years, astronomers have not yet converged on an 



approach for image analysis (Bridle et al. 2009). 



Most image analysis takes one of two approaches: (1) fit an analytical form to the 
light profile such as the Sersic profile to the galaxy distribution or (2) use an orthogonal 
basis function to characterize the ellipticity of an image. While Sersic profiles have proven 
to be remarkably successful at fitting galaxy light profiles outside of the central cores of 



galaxies (Kormendy et al. (2009)), applying them to galaxy images requires non-linear 



fits to the data (Ngan et al. 2009) and computing the effects of seeing is computationally 



demanding. The later approach (e.g., (Bernstein & Jarvis 2002[ )) is computationally simpler 
and mathematically elegant; however, Gaussians are poor approximations to galaxies. 

The goal of this paper is to define a basis function that combines the advantages of both 
approaches. §2 introduces a series of functions, u v K v {u) that have a simple representation 
in Fourier space, so can be easily used in image analysis, and are good approximation to 
Sersic profiles and more importantly to galaxy photometry. §3 generalizes these functions 
to triaxial systems. §4 applies these functions for galaxy photometry and §5 considers the 
effects of seeing and its implications for lensing measurements. 



2. Analytical Functions for Starlight Profiles 

We consider a family of models whose 3-d Fourier profile has the form: 
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where r is the half-light radius, c u is a constant given in table 1 and v > — 1. Projecting 
this to two dimensions, the Fourier transform retains its simple form: 
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Moving to real space, these profiles correspond to analytical 2-dimensional profiles: 
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where 
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and K v {u) is a modified spherical Bessel function of the third kind. The Appendix describes 
some useful properties of these functions. 



For v — j + 1/2, this has a simple form: 
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/i/ 2 («) = exp(-u) 



h/2{u) = -exp(-«)(u+l) 



(5) 

(6) 

(7) 



At small u, /„(«) -»■ l/(2(i/ + 1)) for v > and /„(«) -»■ w 2 7(2(z/ + 1)) for -1 < v < 0. At 
large u, f v (u) — >■ \/tc/2 exp(— ■u)-u _1 / 2+ ' y . 



By going back to Fourier space and integrating by parts, we can analytically evaluate 
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Fig. 1. — This pair of figures shows that a v = —0.6 density profile and a deVa- 
couleur's law profile (m = 4 Sersic profile) have similar forms. The left panel 
compares the surface density profile for the two models. The right panel com- 
pares the integrated surface density profiles(Equation ([8]) 



the integrated light profile of the one component model, 
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Using equation ^8j), we can solve (1 + is)f u +i(c u ) = 1/4 and determine the half-light radius 
for the distribution. Table 1 lists the values of the c u so that r is the half-light radius for 
the light distribution. 



We can also evaluate the 3-dimensional density profile, 
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We can evaluate the potential, 
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Table 1: Half-light Radii 
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by first multiplying by r and differentiating 
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and then integrating again with respect to r to find, 
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where L u is a modified Struve function. 
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For v — 0, this reduces to a simple form:: 
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Since v = is intermediate between an elliptical galaxy profile [y = —0.6) and an 
exponential profile [y = 0.5), it is a simple potential for a plausible generic stellar model. 

We can write the potential in an alternative form by recalling the large x expansion of 



KJx): 
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and one of the definitions of the incomplete Gamma functions, 

poo 
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where T(a,x) is the incomplete Gamma function and Wj = Il|, =1 [(2z/ + l) 2 — (2j — fc) 2 ]. 

There are several interesting extensions to the one component representation: 

(1) Fixed v models An alternative to using the index of the functions (u) as one of 
the parameters in the galaxy fit is to represent the galaxy light profile as a sum of an 
exponential term [y = 0.5) and a very extended profile selected to fit galaxies like M87 
(zy = 0.85): 



E(r/r ) = L\ exp 



+L 2 



-1.6783 
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J J ^0.85(0.35) 

where ro is fit to the scale length of a given galaxy. The coefficient 0.35 was selected so 
that the two profiles span the same ranges as the Sersic functions, so that most galaxy 
profiles can be represented as a sum of the two functions with the same value of ro- Figure 
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Fig. 2. — This figure compare a linear combination of the modified Bessel function 
fit (solid lines), /3/b.5( c o.5 r A'o) + (l — /3)/_ .85(0.35r/r ), to the family of Sersic profiles. 
The solid lines shows the density profiles for fix = 0,0.5 and 1. The dashed lines 
show a set of Sersic profiles with mn = 1,2 and 4. The n = 1 Sersic profile and 
the /5i = 1 profile are identical. 



10 



[2] compares Sersic profiles with m = 1,2, and 4 to surface profiles fit with different linear 



combinations of the two terms in Equation (17). While the two terms are normalized 



to have the same amplitude at r = 7*0, the later term contains more mass for /3 = 0.5. 
Note that this fitting function has three free parameters (Li,L,2,ro), the same number of 
parameters as the Sersic profile (L, n, r ) and the profile introduced earlier in this section 



(L, 1/, ro). The main advantage of Equation (17) is apparent in q4] where we show that we 



can precompute the effects of atmospheric seeing. 
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Fig. 3. — This figure compares Kormendy et al. (2009) compilation of observa- 
tions of M87 (NGC 4468) to their best fit Sersic profile (dashed line with m = 11.8 
and r e = 703.91kpc), a one component profile (solid line with v = —0.85 and 
r o/ci/ = 1944kpc) and a two component profile (cyan solid line with tq/c v = 1944kpc, 
v = -0.85, and r core /r = 0.0043). 
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(2) Core Radii For galaxies with well defined cores but extended light profiles (see 
Figure [3]) , we can generalize the v < profile by adding a second component and an 
additional parameter, r core : 
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This two component form approaches an exponential for r < r core < r . This generalized 
surface brightness distribution also has a simple representation in Fourier space: 



E(k) 
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(19) 



(3) Triaxial Galaxies. These profiles can be generalized to triaxial galaxies by 
introducing a a change of variables to a new set of coordinates, r = Fit and q = F _1 s*, so 
that we can rewrite our profiles, 
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in a new set of coordinates: 



g(Fu) = det (F^ 1 ) / d n s g(F~ l s) exp(ra ■ s) 



(21) 



Applying the transformation, k x — > k x \/l — e and k y — >■ k y ^f{\ +e) the two dimensional 
light distribution in Fourier space becomes, 
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which has a corresponding real space representation: 
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We can expand out the Fourier space profile in a Taylor Series: 



w a T sr r ( 2 + z/ 
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3. The Effect of Seeing and Finite Resolution 

In this section, we use the profiles to explore the effects of atmospheric seeing and 
telescope distortions on determinations of ellipticity. 

The observed image, 1(f), is a convolution of the galaxy's surface brightness profile, 
£(r) with the effects of the atmosphere and the telescope optics:, 

1(f) = d 2 kexp(ik-f)i:(k)R turb (k)R telescope (k) (25) 

where Rturb is the effect of atmospheric turbulence and Rteiescope is the response function of 
the telescope. Using Kolmogorov turbulence theory, 

R(k) = exp(-(kb) 5/3 ) (26) 

where Ofwhm = 2.92076 and b is a parameter that characterizes the correlation length of 
the atmospheric turbulence. Since most astronomical observations are seeing-limited, we 
set Rteiescope = 1 for this section. 

Seeing reduces the signal/noise of the fits of galaxy (and star) profiles to the 
observations. We can estimate the signal to noise by considering fitting a multi-parameter 
model to the data: 

X / 2~[^measuredyi) -1 model \-™-, ^i)\ \ ' ) 

i 

where n® is the noise times unit area, AQ is the area of the pixel, r, is the pixel centroid, 
and Imeasured is the measured light profile and I m odei is a model with a number of parameters 
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A. The noise estimate assumes that we are observing faint galaxies and are dominated by 
atmospheric noise (and/or read noise). 



The simplest case is a point source of intensity A: 



AO, 



X 2 = y_. -^{Imeasured{ri) ~ AZps(rj)] 2 



(28) 



The signal-to- noise for point source detection is: 
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In the continuum limit, we can use Parseval's theorem to evaluate the signal-to- noise: 
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For a circularly symmetric galaxy fit with the profile of section 2, the observed profile 
can again be calculated by the convolution integral: 
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(31) 



where k = kr Q /c u . The effect of atmospheric seeing is a function of b/r or equivalently a 
function of the ratio of the FWHM of the PSF to the effective radius of the galaxy. Figures 
[4] and [5] show that atmospheric seeing "moves" light from the central cusp outwards and 
circularizes the inner portions of the galaxy. These effects are less dramatic for the outer 
isophotes. 
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The signal-to-noise for the galaxy detection is the second derivative of x 2 with respect 

to M: 

/ .<? \ 2 / .<? \ 2 P„ 

(32) 
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where (S/N)ps is the signal-to-noise for the detection of a point source of the same 
magnitude. 
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1 + P 

Note that the effect of galaxy finite size is to degrade the S/N by the ratio of the effective 
areas of a seeing convolved point source to the seeing convolved galaxy. 




0.01 




Fig. 4. — This shows the effects of seeing on the galaxy profile. The four lines in 



the left panel show the amplitude of the monopole term (Eq. 31) for Ofwhm/tq 



0, 0.5, 1.0 and 2.0. The left panel is for an elliptical profile {y = —0.6) and the right 
panel is for an exponential profile {y = 0.5). 
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Seeing has an even more dramatic effect on the ellipticity of the image. The convolution 
of the cos(20) term in the light profile, 
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(34) 



is shown in Figure [5j Seeing makes the central region of the galaxy round and has a reduced 
effect in the outer profile. Because the outer profile is less effected by seeing, it is useful to 
use a functional form that is a good fit to the outer profile for measurements of ellipticity. 
Galaxies are not well-described by Gaussians and the Hermite profile based determination 
of ellipticities are "missing" the information in the outer profiles. 



We can again estimate the signal to noise by differentiating the fit: 



(35) 



Because the cos(20) term is orthogonal to the symmetric term, the error on the second 
term is again just the second derivative with respect to its amplitude: 
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(36) 



where the (1/2) factor comes from angle averaging cos(20 r ). Using Parseval's theorem and 

(37) 



Equations (31- 34) 
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Fig. 5. — This shows the effects of seeing on the galaxy profile. The four lines in 



the left panel show the amplitude of the quadrupole term (Eq. 34) for OpwHM/ r o — 
0, 0.5, 1.0 and 2.0. The left panel is for an elliptical profile {y = —0.6) and the right 
panel is for an exponential profile {y = 0.5) 
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-2 / ^ \ -P2-P0 



This yields: 

^^ l^j P5 (P2 + K)PPS ^^ 

Figure ^ quantifies the additional integration time needed to measure the ellipticity with 
an uncertainty, cr e — 1/5 by plotting the point source signal to noise needed as a function of 
the ratio of the PSF FWHM to the effective radius of the galaxy, 



Figure (P] — fsj) show the effects of atmospheric seeing on lensing measurements. The 
left panels in Figures Q and ^ show that as the seeing degrades, the image gets broader. 
While the right panel shows that as the seeing degrades, the images appear rounder, 
particularly in the inner regions. Note that seeing has less of an effect on the outer 
isophotes. Thus, high S/N observations can recover accurate ellipticity measurements, even 
if 9 fwhm > rQ. However, since the information about the ellipticity is only in the outer 
isophotes, it is important to use an optimized weighting scheme. 



4. Fitting the Profile to Observations 

This section presents an algorithm for fitting light profiles to multi-image stacked data. 
The algorithm utilizes a linearized version of the two component light profile introduced in 



the previous section (see equation (42)). Since the goal of this approach is to have a fast 
linear algorithm that can be applied to large data sets, the profile is expanded in a power 
series in the image size and shape and the convolution of the seeing with each of the terms 
in the series is precomputed. Since we are also expanding the PSF in a series of terms, 
we can simultaneous fit for multi-frame images with little increase in computational cost. 
The final step in the fit is a non-linear step that is not very computationally intensive and 
returns the intensity, size, shape, profile, and orientation of each galaxy image. 
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Fig. 6. — As the seeing degrades, longer integrations are needed to achieve the 
required sensitivity to shear. This figure quantifies this effect by plotting the 
point source S/N need to achieve a statistical error on the shearr of a e ) = 0.2. 
This point source S/N is shown as a function of the ratio of the FWHM of the 
PSF to the effective radius of the galaxy. The solid line is for a v — —0.6 profile, 
which corresponds to an elliptical galaxy. The dashed line is for a v = 0.5 profile, 
which corresponds to an exponential profile. 
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Following the approach outlined in section 2, the galaxy profile is parameterized as 
a five parameter fit: the amplitude of the exponential profile, L 1; the ampltude of the 
elliptical M87-like profile term, L 2 , a flattening term, e, a scale radius r$ and an orientation, 

E(f|Li,L 2 ,r o ,e,0 o ) = Si(f|Li,r o ,e,0 o ) + T 12 {r\L 2 ,ro,e,(j)o) 

t f ( yi + ^os(2(0-0 o ) ^ 
= L 1 a 1 f 1/2 I C1/2 I (41) 

^, f / r nOK Vl + ecos(2(0-0o )' 

+L 2 a 2 f^o.85 0.35— 

where ai = exp(1.6783) = 1-/ f\/ 2 {c\/ 2 ) and a 2 = l//_o.8s(0.35). This fit has a simple 
representation in Fourier space: 

T T^ (42) 

S {1 + ^% 2 [1 " "os(2(0 - O )]} 1+ " 

where z/i = 0.5, v 2 = —0.85, 71 = l/ci/2, and 72 = 1/0.35. Note that we are fitting a single 

value of ro, the effective radius to the prorfie. 

The first step in the analysis is to fit for the position of the galaxy, 9 g and estimate 
its size from measuring its half light radius. By first computing the convolution of a series 
of circularly symmetric galaxy profiles with different characteristic size (e.g., in steps of 
0.1"), the profile fitting calculation becomes a linear problem. Instead of fitting for ro, the 
effective radius, we can fit for 

A S l-g)' ,43) 

where r^ is one of the nearest precomputed fit to the initial list of values . Expanding the 
profile as a Taylor series in (A + (1 — A)ecos(2(0 — 0o)) : 



S(fc|L 1 ,L 2 ,A,,e,0 o ) = Yl 



[1 + fc=7?rg(l - A)(l - ecos(2(0 - O ))] 



l+Vi 



£ £ [A+ (1 - A) £ co S (2(, - ,„))]•< -^f^L (44) 
i T^o [1 + 7^ 2 ?o] 
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where Wj(ui) = U 3 s=1 (i>i + s)/j\. We then rewrite the profile as a sum of terms with linear 
fit coefficents: 

£t(fc|a,) = ^2 a>ijm(Lu L 2, A, e, (p )fx ijm (k, r<) = J^ a^fc, r<) (45) 

jm q 

where the sum over i,j and m is represented as a sum over q for notational simplicity: 

a - m = L^- m {l - A) m e m exp(-2tm<f )0 ), (46) 

and 

g _ j! (7^r i ) 2j w i (z/ i )exp(2im0) 

^' ml ' lj (j-m)!m! [1 + j?k*rlP +j+l " l J 

The next step is to represent the time-varying and spatially varying PSF can be 
expanded as a sum of specified (but not necessarily orthogonal functions), P s (see e.g., 



dJarvis fe Jaih||2004[ )): 

i2 PSF (fc, S , t n ) = ^ £/,(0 fl , t n )P s (k), (48) 

s 

] where t n is the time of each frame, 9 g is the galaxy position, and P s (k) is the Fourier 
transform of the basis function. Typically, U s (8,t p ) will be fit to stars in the image. 

The image profile in a given frame can now be expressed as a sum,: 

I(f, t n \d q ) =y^ y a q y^ y U s (9 g , t n )fjL qa (f, n), (49) 

q s 

over a series of terms that can be precomputed once for a series of effective radii, r,: 

/V(T) = / d 2 kexp(-ik- r)P 8 (k)fi q (k,ri) (50) 

These function, j2ij ms , describe the convolution of a term in the expansion of the galaxy 
light profile with a term in the expansion of the PSF. 

Fitting the convolved linear profile to the observations is now a linear process that fits 
the galaxy light profile coefficients, a^, to all of the data in the stack: 

2 V - ^ {.* measured \J"l j ^n ) 2-^iq a q2-^sUs r g,tn)tJ'qs{ r - l r i))) / K -.\ 

x =2^ ~2 ( 51 ) 
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where the sum over I is a sum over pixels. The best fit is the solution to a linear equation, 

( 52 ) 

where 

(53) 

(54) 

are evaluated as a sum over pixels and frames using precomputed functions. This linear fit 
takes N stack N pix Ng alterms steps to evaluate N kk and N^ alterms steps to invert the matrix. If 
we expand to first order in the ellipticity, then N ga i terms = 6. If we work to fourth order, 
then Ngaiterms = 15. This very rapid process should enable fits to individual observations in 
the stack. 

We then do a nonlinear fit to solve for the amplitude of the terms in the expansion, L\, 
L 2 , e, 0o, and A: 

X 2 ( S |L 1 ,L 2 ,A,e,0 o )= E qq > [af L -af(L l5 L 2 ,A,e,0o)]iV^ 

[oJf L -a£*(£i,L2,A,e^ )] (55) 



where a q is defined in Equation (46). This operation is very quick as it takes only 5N ga i terrns 



steps to evaluate the non-linear fit. 

5. Conclusions and Next Steps 

This note introduces a basis function that may prove useful for the analysis of galaxy 
images, particularly for lensing work. In a subsequent paper, we will test this approach 



against the GREAT08 simulations (Bridle et al. 2009). 
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There are several possible extensions to this parameter fit. We can include priors on 
e and modified the \ 2 (Equation (55) to be a likelihood function and then marginalize 



over the galaxy size and position. We could then also include priors on the ellipticity 



functions (Kitching et al. 2008). Other possible generalization would be to include PSF 



uncertainties and to generalize the expansion to include shapelet terms. 
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A. Appendix: Useful Properties of Bessel Functions and Incomplete Bessel 

Functions of the Third Kind 



Kostroun (1980) provides a useful expression for the numerical evaluation of K v (u): 

exp(— u) 



KJu) = h 



+ \ ^ ex P(~ ucosh(rh)) cosh(z/r/i) 



(Al) 



r=l 

where for values of v of interest to this paper, we set h = 0.5 and find that the sum 
converges to 10~ 10 with less than 7 terms for most values of v and u used in the paper. 



This manuscript was prepared with the AAS TATpX macros v5.2. 
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There are useful relationships between various modified Bessel functions: 

k m+1 dk 



1 



J m [kr) 



/r\ m + u K v (r) 
'I _|_ fc2)i+™+" " \2J r(l + m + z/) 



(A2) 



im(A;r) 



k m+2 dk 



;i + p) 



1+m+i/ 



2/ 



2 






r(l +m + z/) 
^Fr(zy + 3/2) 
+ m + f) 



For z/ < 0, we can use |Abramowitz & Stegun| ([1972) equation 9.6.25, 

exp(— 2 cosh t) sinh v to?t 



/*(*) 



7T 



r(i/ + i)r(i/ + 1/2) 



7T 



exp(— ruj 



-dit 



Thus, 



fu{z)dz 



I> + 1)I> + 1/2) A ( W 2_ i)V2-" 

'vrexp^— x /" exp(— xt>) 



i> + i)i> + 1/2) J (i + v ) („( 2 + 1;)) 1 / 2 - 

Expanding 1/(1 + f )/(2 + ■u) 1 / 2-1 ' in a Taylor series and integrating yields: 

(_i)ir(i/ + |+j) 



(if 



f v (z)dz 



7rexp^— xj 



E 



1 , ^-(1/2-^) 
2ij! 



2"-!/ 2 r(z/ + i)i> + 1/2) ^ (x)V2-^-i 

where Wq(/j,) = 0, W\(n) = /i, w 2 (fi) = 5/x + /i 2 , w 3 (fi) = 32/i + 9/i 2 + /i 3 , and 



(A3) 



(A4) 



(A5) 



(A6) 



w 4 (/x) = 262/x + 83/x 2 + 14// 3 + //. Using |Abramowitz & Stegun| ( |1972[ ) equation 9.6.28 



z dz 



l\ m T(u-m + l) 



r(u + i) 



Jv—m\Z) 



(A7) 



